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Abstract. A rapidly growing amount of evidences, mostly coming from the recent gamma- 
ray observations of Galactic supernova remnants (SNRs), is seriously challenging our un- 
derstanding of how particles are accelerated at fast shocks. The cosmic-ray (CR) spectra 
required to account for the observed phenomenology are in fact as steep as i.e., 
steeper than the test-particle prediction of first-order Fermi acceleration, and significantly 
steeper than what expected in a more refined non-linear theory of diffusive shock acceler- 
ation. By accounting for the dynamical back-reaction of the non-thermal particles, such a 
theory in fact predicts that the more efficient the particle acceleration, the flatter the CR 
spectrum. In this work we put forward a self-consistent scenario in which the account for the 
magnetic field amplification induced by CR streaming produces the conditions for reversing 
such a trend, allowing — at the same time — for rather steep spectra and CR acceleration 
efficiencies (about 20%) consistent with the hypothesis that SNRs are the sources of Galactic 
CRs. In particular, we quantitatively work out the details of instantaneous and cumulative 
CR spectra during the evolution of a typical SNR, also stressing the implications of the 
observed levels of magnetization on both the expected maximum energy and the predicted 
CR acceleration efficiency. The latter naturally turns out to saturate around 10-30%, almost 
independently of the fraction of particles injected into the acceleration process as long as this 
fraction is larger than about 10~^. 
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1 Cosmic-ray spectra: need for steep 

According to the well-known supernova paradigm, the bulk of the cosmic rays (CRs) we 
detect at Earth in a vast range of energies spanning from about a GeV to about lO^GeV is 
of Galactic origin, and supernova remnants (SNRs) are the best candidates as acceleration 
places. 

The reason of the popularity of such a scenario is indeed twofold: on one hand, as 
already put forward many decades ago [1], this class of astrophysical objects can account for 
the required CR energetics and, on the other hand, a very general and efficient mechanism [2] 
turned out be very effective for accelerating particles diffusing around SNR forward shocks 
[3-6]. 

Such a mechanism, which is usually referred to as diffusive shock acceleration (DSA), is 
particularly appealing because it naturally predicts the energy spectrum of the accelerated 
particles to be a power-law oc E~'i, whose spectral index q = does not depend on the 
microphysics of the scattering processes but only on the compression ratio r felt by diffusing 
particles. Moreover, since for any strong (i.e., with sonic Mach number much larger than 1) 
shock r = 4, it is straightforward to predict SNR strong shocks to accelerate particles with 
spectral index q = 2. 

For many years, the most prominent evidence of the non-thermal activity ongoing in 
SNRs has been the detection of synchrotron emission in the radio band due to relativistic 
electrons [e.g. 7], lately corroborated by X-ray observations assessing the extension of the 
electron spectrum up to multi-TeV energies [8]. 

Finally, in the last few years, the present generation of 7-ray instruments has opened an 
additional window on the SNR phenomenology in terms of non-thermal content: the high- 
quality data from GeV-range satellites (AGILE and Fermi), and from TeV-range ground- 
based Cherenkov telescopes (CANGAROO, HESS, VERITAS, MAGIC,...) provided us with 
unprecedented insights into the physical processes involving both accelerated leptons and 
hadrons. We have in fact objects, like SNR RX J1713. 7-3946, whose 7-ray emission is likely 
due to very energetic electrons radiating through inverse-Compton scattering against some 
photon background [9], and objects like Tycho's SNR, whose emission is likely due to the 
decay of neutral pions produced in nuclear collisions between relativistic nuclei (with energies 
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as high as ~500 TeV, at least) and the background plasma [10]. With a terminology that 
has become rather popular, we refer to these scenarios as leptonic or hadronic depending on 
the population of non-thermal particles mainly responsible for the observed 7-ray emission. 

This brand new, and constantly increasing, phenomenology of SNRs can significantly 
probe our comprehension of how DSA works in this class of sources, and even further informa- 
tion is expected to come from the next generation of 7-ray experiments, currently represented 
by the Cherenkov Telescope Array [11]. 

As pointed out for instance in ref. [12], the photon spectral index in most of the 7-ray 
bright SNRs is inferred to be appreciably larger than 2. While in the TeV range this may be 
due to the fact that we are dealing with the cut-off of the parent particle distribution (in both 
scenarios), a rather steep spectrum in the GeV range cannot be ascribed to inverse-Compton 
scattering (which would return a photon spectrum oc E~^'^ for a E~'^ electron one), therefore 
representing a strong signature of hadron acceleration. 

In this paper we do not want to address the question whether most of the SNR 7-ray 
emission is either of hadronic or of leptonic origin, but we want instead to stress a point which, 
in our opinion, has been given too little attention in the literature: 7-ray observations are 
building up more and more evidences that SNRs accelerate particles with spectra different, 
and namely steeper, than what predicted by test-particle DSA at strong shocks. 

In the hadronic scenario, in fact, the 7-ray spectrum has to be parallel to the one 
of parent hadrons. In the leptonic scenario, instead, in order to be consistent with the 
steep spectra (oc E~'^''^-E~'^'^) observed in the GeV band, a contribution from non-thermal 
bremsstrahlung from ~GeV electrons must be added to the inverse-Compton one to fit the 
data [see, e.g., the case of Gas A in ref. 13]. Also for bremsstrahlung emission, however, the 
photon spectrum is parallel to the one of the parent particles, in turn implying that also GeV 
electrons must have a spectrum steeper than E~'^. It is worth recalling that in this region of 
the electron spectrum no cooling via synchrotron emission is effective, therefore protons and 
electrons are expected to show the same spectral index. 

Quite intriguingly, for about 30 years scientists have been working out the extension of 
DSA to the case in which GRs are not simply test-particles but instead participate actively in 
the shock dynamics, carrying sizable fractions of pressure and energy density of the system 
[see, e.g., refs. 14-17, for thorough reviews on the topic]. Such a non-linear theories of 
DSA (usually referred to as NLDSA) invariably predicts the back-reaction of the accelerated 
particles to induce in the upstream the formation of a precursor in which the fluid is slowed 
down because of the pressure in GRs diffusing around the shock. The net result is that 
particles with larger momenta, and in turn larger diffusion lengths, "feel" — on average 
— fluid compressions larger than 4. On the contrary, particles with mildly suprathermal 
momenta only "see" a weaker subshock, with a compression ratio smaller than 4. This spread 
in the fluid compression ratio experienced by GRs leads to rather concave spectra, steeper 
(flatter) than E~'^ at low (high) energies. Moreover, the standard NLDSA theory predicts 
such a concavity to become more and more marked when the acceleration efficiency increases, 
so that fractions of 10-50% of the shock ram pressure channeled into GRs naturally imply 
spectra as flat as E~^-'^-E~^-^ above a few GeV [18]. The steepening induced by the fluid 
precursor can in fact be effective only below this threshold, GeV-nuclei being the carriers 
of most of the GR pressure when the spectrum is steeper than E~'^. Since the photons 
produced in nuclear collisions have a typical energy about 10 times lower than the parent 
hadron, sub-GeV to multi-TeV observations invariably probe energy regions where standard 
NLDSA theory predict spectra systematically flatter than E~'^. 
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No clear-cut evidence of concavity has been found in any of the 7-ray-bright SNRs, but 
this may easily be a consequence of the intrinsic errors in the measurements, which cannot 
be more accurate than a 10-20% in most cases. What was really unexpected, though, is that 
all of the observed spectra [see data collected in ref. 12] are either consistent with, or steeper 
than, the test-particle prediction oc E~'^ and not with NLDSA ones! 

CR spectra as steep as £'~2-3_^-2.4 ^^^^ g^jg^ required when trying to link two distinct but 
actually related aspects of the SNR paradigm, namely the spectra before (at injection) and 
after (at Earth) propagation in the Milky Way. The widely accepted model for CR transport 
in the Galaxy, mainly based on CR isotope and secondary-to-primary compositions, suggests 
the residence time in the Galaxy to scale as E~^, with 6 ~0.3-0.6 [19-21]. The CR flux 
observed at Earth (oc E"'^''^^) has in fact to be proportional to the injection spectrum oc E^'^ 
multiplied by the Galactic residence time oc E~^ , providing the constraint q + 5 ~ 2.75 and 
therefore implying q =2.2-2.4. In addition, the smallest value 5 ~ 0.3 is preferred in order to 
account for the relatively small anisotropy observed in the direction of arrival of CRs above 
~ 1 TeV [see, e.g., 22]. Again, a hint that the spectra of CRs accelerated in SNRs have to 
be non-negligibly steeper than the test-particle case, and significantly steeper than standard 
NLDSA predictions. 

In this respect, it is however important to remember that the spectrum injected into 
the Galaxy during the whole SNR lifetime has not to be compared one-to-one with 7-ray 
observations, since such a spectrum is expected to be a rather complex superposition of 
several time-dependent contributions, as illustrated in ref. [23]. Despite of our still incomplete 
understanding of how accelerated particles leave their source and become CRs, physically 
motivated calculations [24-26] indicate that the cumulative spectrum (i.e., integrated over 
the SNR lifetime) turns out to be only ~ 0.05-0.1 steeper in slope than the instantaneous 
spectrum achieved in the early Sedov stage. 

We can safely summarize the points above by saying that the observational evidences 
are better recovered if the spectra of accelerated particles at SNR shocks were steeper than 
the test-particle prediction oc E~'^ . How this is actually achieved is a fundamental and chal- 
lenging question whose answer cannot be found in standard NLDSA theories, since they 
invariably predict that the more efficient the acceleration, the flatter the spectrum of accel- 
erated particles. 

Therefore, we are left with an open dilemma: are Galactic CRs efficiently accelerated 
in SNRs? And in this case, what is the physical mechanism able to reverse the standard 
prediction of the NLDSA theory? Understanding why observed 7-ray spectra from SNRs are 
so steep is indeed a fundamental tile in the comprehension of particle acceleration at shocks, 
and not an almost negligible correction to an otherwise solid physical theory. 

In this paper we investigate one possible mechanism, put forward already at the dawn 
of DSA theories [5], which, in a revised form accounting also for the more recent evidences of 
magnetic field amplification in young SNRs, has been proposed as a viable scenario for accel- 
erating particles with rather steep spectra consistent with 7-ray observations [12]. The main 
idea here is that, since efficient CR acceleration may induce a very efficient amplification of 
the magnetic field via some plasma instabilities, the magnetic structures acting as scattering 
centers for the CR diffusion may achieve a non-negligible velocity with respect to the back- 
ground fluid [also see refs. 27-30, for other recent implementation of the same effect]. This 
phenomenon may significantly alter the actual compression ratio felt by accelerated particles 
and, in turn, their spectrum and eventually the global shock dynamics. 

The original contribution of this work is a quantitative analysis of such an effect, as- 



-3- 



sessing — as self-consistently as possible — the spectral steepening as a function of the CR 
acceleration efficiency, here regulated by the number of particle extracted by the thermal 
bath and injected in the acceleration process. 

The plan of the paper is as follows: in section 2 the details of the model are illustrated. 
The solutions of the equation describing the CR transport and the saturation of the magnetic 
field amplification are worked out, along with their non- linear interplay with the SNR hydro- 
dynamics. In section 3 we discuss our main results, namely the fact that including the finite 
velocity of the scattering centers in the amplified magnetic field it is possible to reproduce 
many observational features, like the levels of magnetic field amplification, the steepness of 
the CR spectra and the achievement of maximum energies consistent with the knee. We also 
outline the relationships between instantaneous and cumulative spectra in different stages 
of the SNR evolution and, quite interestingly, point out for the first time how our findings 
does not depend in a critical way on the fraction of particle injected, as long as this is larger 
than ~ 10"^. After some comments about strong and weak points of the present model, we 
conclude in section 4. 

2 The kinetic model for particle acceleration 

The final goal of any approach to the problem of NLDSA is to solve self-consistently the 
equations for mass, momentum and energy conservation, along with a description of the non- 
thermal particles and, possibly, their generation of magnetic waves via plasma instabilities. 

Apart from computationally very expensive particle-in-cell (PIC) simulations of colli- 
sionless shocks, in which the interplay between particles and fields is calculated from first 
principles, the most common way to describe the transport of relativistic particles is either 
to prescribe a scattering law [Monte Carlo approaches, see for instance refs. 16, 31, 32] or to 
solve a Vlasov-like diffusion-convection equation for the isotropic part of the CR distribution 
function [33]. More precisely, the diffusion-convection equation has been solved numerically 
in its time-dependent formulation [34-39] or semi-analytically, i.e., by integrating the quasi- 
stationary equation in order to have an implicit analytical expression for the CR distribution 
function [40-45]. Nevertheless, all of the different approaches involving a model for the par- 
ticle transport lead to very consistent results for non-relativistic shocks, as demonstrated in 
ref. [18], and at the moment represent our best description of particle acceleration in SNRs. 

In this work we apply the semi-analytical formalism for NLDSA put forward in refs. [44, 
46] , in its extended version including particle escape from an upstream boundary [45] and the 
dynamical feedback of self-generated magnetic fields [28, 47]. Again, we would like to stress 
that the limitations and the uncertainties in applying this NLDSA model to SNR shocks 
are shared by virtually any non-PIC approach: the main advantage of a semi-analytical 
formalism is to be very quick and versatile (a run takes several seconds on a standard laptop), 
and therefore well-suited for investigating a large parameter space. 

The stationary mass and momentum conservation equation for a plane, non-relativistic 
shock simply read 



where p, u and p represent the plasma density, velocity (in the shock reference frame) and 
pressure, while the subscripts cr and B label the pressure in the shape of non-thermal par- 
ticles and magnetic fields. Throughout the paper the subscripts 0, 1 and 2 refer to physical 
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quantities measured at upstream infinity, immediately upstream of tlie shock and down- 
stream, respectively. The shock is stationary at a; = and, by definition, uq = Vsh, i-e., the 
shock velocity in the Earth reference frame. 

Global conservation of energy is then granted by taking into account the proper jump 
conditions at the shock [47] and by solving the space-dependent equations for the energy 
transport of thermal gas, magnetic waves and CRs. In particular, we consider the upstream 
thermal plasma to be adiabatic, namely 

pixp pI 3 

[a generalization to non-adiabatic heating is straightforward, see ref. 28] and we model the 
growth and the advection of the Alfven modes generated by accelerated particles via stream- 
ing instability by solving the transport equation for the magnetic turbulence [e.g., 65] ac- 
cording to the fluid formalism in ref. [28]: 

dpB{x) .dpcr{x) du{x) 

This stationary equation simply describes the advection of the magnetic pressure (left-hand 
side) as produced by the CR pressure gradient (first term on the right-hand side), also ac- 
counting for adiabatic compression in the precursor (second term on the same side). Here we 
consider Alfven modes generated by CR streaming instability and thus propagating against 
the fluid, therefore u{x) = u{x) + VAi^) is the fluid velocity in the wave reference frame, 
where particle scattering is elastic. 

We will comment the generality of our choice more widely in the following sections, but 
it is worth stressing that eq. 2.4, in its simplicity, still captures a basic feature of the interplay 
between CRs and magnetic field: the CR diffusion velocity ~ ~ « is typically larger 

than the Alfven velocity, therefore the wave-particle coupling tends to restore an equilibrium 
between the two by making va larger and the gradient in pcr smaller, i.e., by amplifying the 
magnetic field and smoothing the precursor. 

Finally, we include the diffusion-convection equation for the isotropic part of the CR 
distribution function f{x,p) [see, e.g., ref. 48]: 



dfix,p) d 

— ^ = IT 

ox ox 



where Q{x,p) accounts for particle injection and 



D{x,p) = ^rUx,p) (2.6) 

is the Bohm-like parallel diffusion coefficient for a particle with velocity v{p) and Larmor 
radius rL{x,p) = in the local, amplified magnetic field B{x). 

Following the implementation in ref. [49], we assume that, immediately behind the shock, 
all the particles in the Maxwellian tail with a momentum larger than a critical pinj have a 
Larmor radius large enough to be able to cross the shock and return upstream [thermal 
leakage, see e.g., 37]. Such an injection momentum is parametrized as a multiple of the 
downstream thermal momentum Pth2^ namely 



Pinj = iinjPth,2; Pth,2 = V^rripkBT^, (2.7) 
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where is the proton mass, the Boltzmann constant and T2 the downstream tempera- 
ture. We thus have 



where 



Qix,p) = r] 



V 



PlUl 



-S{p - Pinj)5{x) 



(2. 



(2.9) 



represents the fraction of the particles crossing the shock injected in the acceleration process 
and 6{x) accounts for injection to occur at the shock position. Since the additional pressures 
in CRs and magnetic fields induce a velocity gradient in the upstream (precursor), it is 
convenient to introduce two distinct compression ratios as felt by scattered particles, the 
subshock and the total ones: 



R-sub — 



Ui + VA,1 _ 



R 



■tot 



Up + VA,0 
U2 



(2.10) 



These compression ratios differ from the fluid ones because of the finite velocity of the scat- 
tering centers va, which vanishes downstream because of wave isotropisation {va,2 ~ 0). For 
typical interstellar values also va,o/uo *C 1, but when efficient magnetic field amplification 
occurs in the precursor va may become a non-negligible fraction of ui [also see the discus- 
sion in ref. 28]. Moreover, when self-generated by CRs, waves tend to propagate in such a 
way to smooth the CR pressure gradient out, i.e., against the fluid: very generally, we have 
u = u + VA ^ u, so that the compression ratios in eq. 2.10 turn out to be smaller than their 
fluid counterparts (see section 3.3 for additional comments). 

The solution of eq. (2.5) with the spatial boundary condition /(xo,p) = 0, which mimics 
the presence of an upstream free escape boundary placed at x = xq, can be written as [45]: 



f{x,p) = /2(p)exp 



df 

^esc{p) = -D{xo,p) — 



W{x,p) 



dx' 



Uq 



D{x>,p) 



exp 



D{x',p) 
_uoMpI 
Woip) 




1 



W{x,p) 
Woip) 



dx" 



u[x 



D{x",p) 



(2.11) 
(2.12) 

(2.13) 



Here $esc(p) represents the flux of particles escaping the system at xq [50] and the 
distribution function at the shock f2{p) = /i(p) reads: 



where 



and 
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/2(P) 



«1 
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A-Kp, 



pip) r dp' 



Upip') + 



dx ( duix) 



U(j I dx 



■ exp 



Dix',p) 




Qpip) 



3R 



tot 



(2.14) 



(2.15) 



(2.16) 



RtotUpip) - 1 

correspond to the average fluid velocity and spectral slope relevant for a particle with mo- 
mentum p. In the test-particle limit, in fact, = and in turn Up = 1 and qp = 
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Even if the calculations are led in the more natural framework of the momentum space 
for the CR distribution function, f{p), in the rest of the paper we will use the energy distri- 
bution function 4>{E) = Anp^ f{p)^ to provide a clearer comparison with observations. 

2.1 Magnetic field amplification 

When particle acceleration is efficient, magnetic field amplification due to CR-induced in- 
stabilities may eventually lead to 6B/B ^ 1 non only in the downstream but even in the 
upstream: such a scenario is motivated, for instance, by the lack of detection of X-ray emis- 
sion from the precursor of SN 1006 [51]. This evidence is to some extent complementary to 
the well-know evidences of large magnetic fields in the downstream, and strongly suggests 
that magnetic field amplification is induced by CRs rather than by some process occurring 
at or behind the shock only. 

Even if the resonant Alfven modes produced by streaming instability [5] are purely 
transverse, the presence of other non-resonant modes and the likely onset of a certain degree 
of turbulence might lead to an — at least partial — isotropisation of the amplified magnetic 
field, in turn affecting the component of the field parallel to the shock normal and eventually 
increasing the effective Alfven velocity, as demonstrated by the authors in ref. [52]. In order 
to heuristically take this effect into account, we assume that in eq. (2.4) the relevant Alfven 
velocity is the one in the local amplified field, namely 



va{x) 



B{x) 
\/47r/j(x) ' 



Bix) 



V8 



ttpb{x) 



(2.17) 



Normalizing the velocities to the shock velocity uq and the pressures to the ram pressure 
PqUq (indicated with capital letters) we can write eq. (2.4) as 



2U{x) 



dx 



Va{x 



dPcrjx) 

dx 



SPb(x 



dU{x) 

I 

dx 



(2.18) 



Since Pb oc eqs. 2.4 and 2.18 are accurate at the second order in Va/U^ therefore the 
correction due to the finite velocity of the scattering centers in the advective and compres- 
sional terms can be neglected and we can approximate, in this equation, U{x) as U{x). 

Moreover, since SNR shock typically show both sonic and Alfvenic Mach numbers much 
larger than one, at the order ^> 1, ^ 1, from eq. (2.2) we have Per — 1 — U{x) and 
therefore eq. (2.18) becomes: 



dPBjx) 
dx 



Pb{x) , 3Pi?(x) 



+ 



2U{x) 2 U{x) 



dU{x) 
dx 



(2.19) 



Introducing ^'(x) = \/ 2u{t) the simple equation 



d^{x) _ 5^'(x) + 1 dU{x) 



dx 



4C/(x) 



dx 



(2.20) 



whose solution with the boundary condition ^(xq) = Pb{xq)/2 ~ is 

1 



1 - U{x) 



-5/4 



(2.21) 
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thus finally leading to: 

This equation, besides taking into account the adiabatic compression of the modes 
in the precursor due to the factor C/(a;)~^/^, connects the local amount of magnetic field 
with the local modification induced by CRs, and therefore naturally predicts the maximum 
amplification to be achieved immediately upstream of the subshock. If this were rigorously the 
case, it would be unlikely to have non-linear magnetic structures with large phase velocities 
far in the precursor, where the high-energy particles diffuse, and therefore the standard 
NLDSA prediction for the total compression ratio felt by the highest energy particles would 
be recovered. 

On the other hand, there are many effects that are not included in eq. (2.18) which may 
play a relevant role, like for instance the excitation of Bell's non-resonant modes [53, 54], 
or the development of large scales structures due to fire-hose instability [55, 56], which may 
be most effective where the CR distribution function is most anisotropic, i.e., far in the 
precursor. There are also phenomenological reasons for the plausibility of such a scenario: 
in Tycho's SNRs there are evidences for macroscopical magnetic structures on the scales of 
the Larmor radius of ~ lO^GeV particles [57]. 

In the absence of a consistent kinetic theory for the CR-magnetic field interplay in the 
precursor, we will assume in the rest of the paper a constant magnetic field in the upstream, 
whose strength is given by the saturation at the shock from eq. (2.22), a choice that proved 
itself to be adequate in particular for the account of Tycho's multi-wavelength emission and 
hydrodynamics [10]. We will further comment on this important point in section 3.3. 

Eq. (2.22), coupled with eq. (2.3) for the background plasma and the solution of eq. (2.5) 
for the CRs, allows for global energy conservation. The connection between subshock and 
total compression ratios for the fiuid can be finally worked out from the solution of the 
Rankine-Hugoniot at the subshock also including the dynamical effect of the amplified mag- 
netic field according to equation 16 in ref. [28]. 

A full solution of the problem of NLDSA can be obtained by solving iteratively, at any 
given time, the system of equations including also eqs. (2.1) and (2.2) through the procedure 
put forward in ref. [45]. 

2.2 SNR evolution 

The SNR evolution is modeled following the analytical recipe given in ref. [58] , and in this pa- 
per we consider the propagation of the forward shock in a homogeneous medium with particle 
density no = Po/mp = O.Olcm"^ and temperature Tq = lO^K except when otherwise stated. 
The total kinetic energy and mass in the ejecta is fixed in Esn = lO^^erg and Mgj = 1.4M0, 
with Mq the solar mass. The generalization to a more complex circumstellar environment 
including winds and bubbles produced in the pre-SN stages is rather straightforward [12], 
but it is omitted here for clarity's sake. 

With the environmental parameters chosen the forward shock dynamics is described by: 

Rsh{t)c^U.1t'/''pc; Vsh(.r)c^4U0T~'/''— (2.23) 

s 

during the ejecta-dominated stage (r = t/Tsr < 1), with Tst — 1900yr, and by 

Rsh{r) ^ 16.2 (r - 0.3)^/^ pc; V.hit) ^ 3330 (r - 0.3)"^/^ — (2.24) 

s 
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Figure 1. Time evolution of several physical quantities: shock radius, velocity and sonic Mach 
number {Rsh,Vsh and M^^q), downstream magnetic field (B2) and CR acceleration efficiency at the 
shock (Per) for ^inj = 3, no = O.Olcm"'^ and Tq ~ fO^K. With these environmental parameters the 
Sedov-Taylor stage begins around 2000 yr and the sonic Mach number drops below f around lO^yr. 



during the Sedov-Taylor stage (r = t/Tsr > 1), as in figure 1. The quasi-stationary solution 
is calculated at many times from r = 0.03 to the end of the Sedov stage, namely r ~ 60 
according to eq. 13 in ref. [58]. The spectra of advected and escaped particles are obtained 
by convoluting the instantaneous spectra weighted with the shell volumes competing to each 
time step as described in ref. [23] . Also adiabatic losses due to the shell expansion are taken 
into account: more precisely, at any time t > to the energy E{t) of a particle with energy Eq 
advected downstream at time to reads [section 3 of 23] 



E{t) = Eo 



Vshjt) 
Vsh{to) 



2 



37 

(2.25) 



with 4/3 < 7 < 5/3 and a rather week dependence on 7. 

Within this framework, we can study the SNR evolution in terms of particle acceleration 
across the transition between the ejecta-dominated and the Sedov-Taylor stages, when the 
CR maximum energy is achieved [59]. 

3 A modern view of CR-modified shocks 

In this section we want to quantitatively investigate the non-linear response of the system to 
the presence of accelerated particles and magnetic fields. In order to do this, we also need to 
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specify the fraction of particles injected in the acceleration process, rj. This quantity has not 
been worked out from first principles, yet, and the only insights come from PIC simulations 
of collisionless shocks in given regimes [e.g., 60, 61]. Since in our framework, as in any other 
approach aimed to account for the temporal and spatial scales relevant to SNRs, we do not 
have any constraint on the microphysics that regulates particle injection, we allow for the 
variation of $,inj, actually the only free parameter in our model. 

In the context of the thermal leakage model rj increases for smaller ^inj, i.e., when the 
minimum momentum required for a particle to cross the shock from downstream is closer 
to the thermal momentum (eq. 2.9). More precisely, in this paper we will span the range 
r/ w 5 X 10~^ — 5 X 10"^ by changing between 4.5 and 2.5. 

3.1 Instantaneous and cumulative spectra 

The evolution of CR spectra is illustrated in figure 2, where the instantaneous spectra of 
both advected (thick) and escaping particles (thin lines) are shown, xq = O.bRsh is chosen 
throughout the paper. The two earliest time-steps depicted fall in the ejecta-dominated stage, 
while the others describe the early, intermediate and late Sedov-Taylor stages. The top and 
bottom panels correspond to two injection efficiencies: (,inj =3 and 4, i.e., rj ~ 2.5 x 10~^ 
and 7] = 4.1 x 10~^, respectively. 

There are many points worth noticing in discussing figure 2. 

• As expected, the energy channeled into CRs scales with the number of injected particles: 
the smaller (larger rj), the larger the energy in non-thermal particles, as can be 
inferred by comparing the normalization of the curves in the two panels, at any time. 

• As a consequence of the larger CR pressure, the case = 3 provides a significantly 
larger magnetic field amplification, and therefore the achievement of a significantly 
larger maximum energy in the early evolutive stages. More precisely, in the efficient 
case (top panel) the spectra turn out to be cut-off around 3 — 5 x 10^ GeV, in nice 
agreement with the position of the knee of the diffuse spectrum of Galactic CRs. 

• For a given injection efficiency, the CR spectra are almost constant during the ejecta- 
dominated stage and extend up to the maximum energy achievable during the SNR life. 
Nevertheless, when the shock begins to slow down appreciably because of the inertia 
of the swept-up mass, the instantaneous maximum energy starts to decrease. The 
other main effect of the velocity drop is the reduction of the Alfvenic Mach number, 
which leads to steeper and steeper spectra. The CR acceleration efficiency, however, 
remains almost constant as long as the shock is still strong (M^^o ~ 10 up to about 
lO^yr, as depicted in figure 1). Such a trend at intermediate ages is non-trivial, and 
is a consequence of the dependence on the fluid velocity of the streaming instability 
saturation (see eq. 2.22). In other words, the effective Alfven velocity induced by 
magnetic field amplification decreases more slowly than the shock velocity, therefore 
leading to a smaller Alfvenic Mach number and, in turn, to smaller compression ratios 
felt by the accelerated particles. The steepening of the spectra and the decrease of 
Emax turn out to be actually compensated by an increase of the number of particles 
at lower energies, say around 1-10 GeV. At the end of the Sedov-Taylor stage also 
the sonic Mach number drops. Per decreases almost linearly with time and the spectra 
become even steeper. 
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Figure 2. Instantaneous CR spectra for different SNR ages as in the legend. Thick (thin) hnes show 
the spectra of particles advocated downstream (escaping from the upstream free escape boundary). 
The panels correspond to two different injection parameters: = 3 (top) and ^inj ~ 4 (bottom 
panel), corresponding to 77 ~ 2.5 x 10^'^ and 77 = 4.1 x 10~^, respectively. 
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Figure 3. Different curves represent tlie downstream magnetic field produced by CR- induced stream- 
ing instability in young SNRs for different injection efficiencies and circumstellar particle densities, as 
in the legend. The data point for historical young SNRs are worked out from the best fitting of X-ray 
bright rims as due to synchrotron losses, with a fiducial error of 20% [see refs. 10, 47, 62-64]. 



• Despite of what predicted by standard NLDSA theories, the more efficient the CR 
acceleration, the steeper the spectra of the accelerated particles (compare the cases 
Cmj = 3 and 4 in figure 2). This effect is a consequence of our assumption of calculating 
the velocity of the scattering centers in the amplified magnetic field rather than in the 
background one, and can be understood in terms of the effectiveness of the magnetic 
field amplification. We will comment more widely on this point in section 3.2. 

The present findings therefore rely on the magnetic fields to be produced by CR-induced 
plasma instabilities and more precisely on the presence of modes propagating upstream with 
respect to the fluid with an effective phase velocity close to the Alfven speed in such enhanced 
fields. 

It is important to stress that, for the more efficient case £,inj = 3, the self-generated 
magnetic fields that produce such a steepening are consistent with the ones inferred by X- 
ray observations. In figure 3 the predicted magnetic fields for different injection efficiencies 
and circumstellar densities are in fact compared with the downstream fields in very young 
Galactic SNRs, as inferred by interpreting the width of their X-ray-bright rims as due to 
synchrotron losses [10, 47, 62-64]. 

It is worth recalling that the NLDSA theory is actually able to predict the CR accel- 
eration efficiency and the pressure in magnetic fields as fractions of the bulk pressure poV^j^, 
therefore the absolute value of B2 is rather dependent on our knowledge of the SNR hydro- 
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dynamics and of the circumstellar density. This effect is shown in figure 3, where a larger 
no naturahy leads to a larger B2, at fixed ^inj, while at fixed particle density the magnetic 
field increases with the injection efficiency. This degeneracy can be broken only if detailed 
information about no and shock velocity is available through a model of the SNR evolution, 
which in turn requires knowing its age and distance, and/or through the detection of some 
hadronic emission, which is related to the absolute value of po- SNRs for which all of these 
constraints are available are very rare: at the moment Tycho's SNR represents the best case 
in this respect [10]. 

Some of the effects outlined above are also visible when the cumulative (as opposed to 
the instantaneous) CR distributions are plotted, as in figure 4. The thick curves correspond 
to the total number of CRs advected in the downstream up to the age in the legend, also 
accounting for the adiabatic losses due to the expansion of the shell (eq. 2.25). 

The thin, bell-shaped curves in figure 4 show instead the cumulative spectra of the 
particles which have escaped the SNR from the upstream boundary. It is easy to see that 
escape always occurs at the highest energies, and that the total escape flux grows rapidly 
during the early stages, reaching a saturation around ~ 10 kyr for the parameters chosen 
here. The escape flux is expected to be larger when the CR spectrum is flatter [50], and in 
fact such a trend is consistently recovered in figure 4, where the normalization of the escape 
flux is larger in the bottom panel (case ^inj = 4), despite of the smaller CR acceleration 
efficiency. 

In both cases, the total SNR content in non-thermal particles reaches a sort of saturation 
after a few thousands years: the contribution of particles accelerated at later and later stages 
(with steeper spectra) becomes less and less important and is more and more relegated to 
the lower energies. At large energies, the contribution is not enough even to balance the 
adiabatic losses due to the shell expansion: the net effect is a mild steepening of the overall 
spectrum, visible for both the efficiencies in figure 4. Eventually, adiabatic losses become 
faster than the particle supply at any energy and the total spectrum normalization start to 
decrease (compare the ~ 10 and ~ 30 kyr curves). 

Quite interestingly, since the CR acceleration efficiency turns out to be almost constant 
during the early stages up to about lO^yr, some of the predictions put forward in ref. [23] 
are recovered also in the present framework. More precisely, the convolution over time of the 
escaped particles turns out to be a power-law spanning the range of variation of Emax] it 
especially accounts for the CRs the SNR can inject in the Galaxy above ~ lO^GeV. Particles 
with the same energy which have been advected downstream suffer relevant adiabatic losses, 
instead: without escape from upstream, SNRs cannot account for the knee observed in the 
diffuse spectrum of Galactic CRs. 

It is very important to stress that, even if the instantaneous spectrum may become 
steeper than during the Sedov stage, the total CR content of the SNR, which is built 
up mainly by particles accelerated when the shock is faster, is not expected to immediately 
reflect such a change. In other words, for middle-aged SNRs, a meaningful interpretation 
of the 7-ray data requires a time-dependent treatment, in that simple estimates of the SNR 
content in CRs according to the instantaneous shock velocity and acceleration efficiency may 
easily lead to both an underestimate of the emission and an overestimate of the spectral 
slope. 

For these reasons, observing spectra significantly steeper than E~'^ in the early Se- 
dov stage implies that consistently steep spectra must have been produced since the ejecta- 
dominated stage: this fact puts an important constraint on the details of how NLDSA works 
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Figure 4. Cumulative (i.e., total) CR spectra at different SNR ages as in the legend, divided in 
advected (thick) and escaped particles (thin lines), for S^inj = 3 and 4 as in figure 2. 
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Figure 5. Top panel: downstream magnetic field {B2, solid red line) and CR pressure at the shock 
{Per, dashed blue line) as a function of the fraction of injected particles 77 (see also the corresponding 
values of (.mj, green triangles). Bottom panel: spectral slope of the particle distribution q, always as 
a function of 77. All the quantities are calculated at the same SNR age, 3500 yr, when there is no 
much difference between the slopes of the instantaneous and cumulative CR distributions. 



in fast shocks. Also for the cumulative CR spectra the same consideration made for the 
instantaneous spectra holds: the more efficient the CR acceleration, the steeper the resulting 
spectra and the larger the maximum energy achievable by accelerated particles (compare the 
case ^inj = 3 and S^mj = 4 in the panels of figure 4) . 

3.2 Injection efficiency 

We want to investigate now in more detail the role of the injection parameter in the deter- 
mination of the expected CR acceleration efficiency. In particular, it is interesting to check 
whether the steepening of the spectra which has been observed for low (large r]) may 
lead to a saturation of the pressure in CRs at the shock. The opposite limit, the one with 
large ^mj (small r/) is somehow less interesting, since it should simply recover the test particle 
prediction and therefore lead to very inefficient magnetic field amplification. 

Let us consider a fixed SNR age of 3500 yr (early Sedov-Taylor stage), and let 
vary between 2.5 and 4.5, corresponding to the injection of a fraction between 5 x 10~^ and 
5 X 10~^ of the particles crossing the shock at any time. 

In the top panel of figure 5 the fraction of the bulk pressure in CRs, Per = Per, 2, and the 
downstream magnetic field, B2, are shown as a function of the fraction of injected particles. 
The value of ^inj that needs to be chosen to inject the correspondent fraction rj is showed in 
the top axis. 
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For rj smaller than ~ 10~^ the acceleration efficiency is less than a few %, and in turn 
magnetic field amplification is very inefficient: downstream we are basically left with the 
component parallel to the shock normal only: Bq ~ From r] ~ 10~^ to r/ ~ 10~^ Per 

increases almost linearly (up to almost 10%) and so does the magnetic field. 

The most interesting effect, though, is that for r] larger than 10~^ the pressure in 
CRs remains almost constant, between 10 and 30%. The variation of almost three orders 
of magnitude in r/, in fact, reflects in a change of less than a factor three in Per and i?2- 
The very reason for this effect is accounted for in the bottom panel of figure 5, where the 
spectral slope of the CR distribution at the shock q = ——^^^^ is shown. When the shock 
modification reaches about 10%, the feedback of the finite velocity of the scattering centers 
kicks in, leading to a steepening of the spectrum with respect to the oc E""^ prediction of 
the test-particle limit. The less and less energy that goes into the high-energy part of the 
CR spectrum (which extends to higher energies, though) is rather fairly balanced by the 
more and more energy channeled into the lower energy particles, in such a fashion that the 
overall pressure and density energy in CRs turn out to be rather insensitive to the injection 
details. Nevertheless, the spectral slope of the CR distribution increases with the increase 
of the injection efficiency up to g ~ 2.6, spanning values consistent with most of the 7-ray 
bright SNRs [12]. 

Such an effect gives a new flavor to the notion of a CR-modified shock: an efficient CR 
acceleration does not lead to an arbitrary increase of the pressure in CRs and therefore to 
spectra asymptotically as flat as oc E'~^-^ [see, e.g., 46, and references therein], but rather to 
a significant steepening of the particle spectra with a moderate fraction (10-30%) of the bulk 
energy channeled into non-thermal particles. This is probably the most interesting original 
finding of the present work, and represents an important step in order to consistently account 
for the new wealth of information coming from 7-ray observations and the long-lasting SNR 
paradigm for the origin of Galactic CRs. 

We want also to stress that other mechanisms as turbulent heating [see e.g., 65, 66] 
and dynamical magnetic feedback [47] have been proposed in order to avoid the arising of 
very modified shocks, and consequently the appearing of very concave spectra as flat as 
j^;-i-2 _ £;-i-5 |;]^g highest energies, which would be at odds not only with 7-ray but also 
with radio and X-ray observations of SNRs. The first mechanism is however unlikely to 
be very effective, in that it would require a severe dissipation of the magnetic turbulence, 
contrary to the several evidences of amplified magnetic fields found in young SNRs. The 
latter, instead, is expected to be rather relevant in the current framework (and it is in fact 
included in the present calculation) and may effectively contribute to preventing an excessive 
shock modification [28]. 

However, neither the first nor the latter piece of Physics, alone or together, can account 
for the required steepening of the CR spectra in that, if one does not allow for a distinct 
phenomenology of the compression ratios felt by the fluid and the CRs, the predicted spectral 
slope cannot in any case be larger than 2 for strong shocks at energies above a few GeV, i.e., 
in the region of interest for 7-ray astronomy. 

3.3 Some additional comments 

A natural question about the flndings illustrated in the previous section is how strong is 
the dependence of the steepening on the (still) not-well understood nature of the magnetic 
turbulence. Quantitatively speaking, for the present mechanism to be effective, the magnetic 
irregularities need to have a phase velocity not-negligible with respect to the fluid one, on a 
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length-scale comparable with the diffusion length of the most energetic particles. Va/U ~ 
10 — 15%, corresponding to effective Ma ^ 6 — 8, is enough to return a spectral slope of 
2.2-2.35. The sign of the correction in the compression ratios (eq. 2.10) comes for free when 
the turbulence is generated by the CR streaming, since upstream waves should propagate 
against the CR gradient, i.e., against the fluid, while downstream it seems unlikely to have 
the propagation of waves with any given helicity since the medium is hot and turbulent. 

Another fundamental property required for the overall steepening of the CR spectra 
is that also particles diffusing far in the precursor must scatter against magnetic structures 
moving with a phase velocity comparable to the ones of modes propagating closer to the 
subshock. It is in fact possible to show that by assuming a local magnetic field amplification 
as given by eq. (2.22) as a function of x, which would imply less and less turbulence closer 
and closer to the free-escape boundary, the CR spectrum remains as flat as E~'^ close to the 
maximum energy. In other words, either energetic particles diffusing up to the upstream edge 
of the precursor decouple from the fluid and escape the system (in a fashion more complex 
than the one described by the transport equation 2.5), or there needs to be enough ongoing 
magnetic field amplification even rather far upstream, on the scales relevant to the highest 
energy particles. 

Whether such a magnetic configuration is really achieved in SNR shock cannot be 
addressed with analytical techniques, since it requires the investigation of very complex 
particle-wave interactions, well beyond the quasi-linear limit. The most useful insights into 
these details come instead from PIC simulations, in which the problem is tackled from basic 
principles [see e.g., 52, 67, 68]. 

Bell modes are expected to provide an additional channel for amplify the magnetic 
field, especially close to the escape boundary where the CR distribution function is rather 
anisotropic. As pointed out in refs. [67, 68], while in the linear regime of the instability 
modes are almost purely growing, the saturation of Bell modes in the non-linear regime 
{5B/Bq ~ 10) may occur exactly because of the local Alfven velocity to become comparable 
with the drift velocity of the CR current. 

An additional promising class of Alfvenic modes is represented by the left-hand, circu- 
larly polarized ion-whistler waves found in hybrid (kinetic protons/fluid electrons) simulations 
described in ref. [61]. In their work, in fact, the authors find that, in the upstream parallel 
shocks with Ma ~ 10, modes with phase velocity larger than the Alfven velocity can be 
efficiently excited. 

Finally, the possibility of producing large-wavelength modes through the fire-hose mech- 
anism has been put forward by some authors when investigating the development of the beam 
instability [55, 56]. 

Even if a lot of work has still to be done in order to understand which instability 
is the most relevant for SNR shocks, these scenarios are quite consistent with our working 
hypothesis that, when some channel for CR-induced magnetic field amplification is active and 
leads to sizable perturbations with 5B/B ^ 1, the wave-particle coupling tends to reduce 
the gap between the velocities of the magnetic structures and of the CR drift: the net effect 
is therefore to force the magnetic structure to acquire a non-negligible velocity with respect 
to the fluid, as heuristically implemented in section 2. 

From the point of view of the phenomenology of 7-ray bright SNRs and their steep 
spectra, it is important to recall that there is another mechanism which may be relevant 
to reconcile theory and observations. As pointed out in ref. [69], when a shock propagates 
into a partially ionized medium, charge-exchange processes between neutral and protons in 
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the downstream lead to the formation of a neutral return flux which eventually conveys 
some of the downstream energy back in the upstream, producing a neutral-induced precursor 
on a scale comparable with the ionization/charge-exchange mean free path (~ 10^'' cm for 
typical SNR parameters). Such a neutral-induced precursor, which may reduce the sonic 
Mach number down to 2 — 3 before the subshock, develops on a scale independent of the 
CR diffusion length, therefore can in principle induce a sizable modification in the velocity 
profile and in turn on the expected compression ratios, possibly leading to steep CR spectra 
up to the highest energies. 

The only limitation of this unavoidable mechanism is that it is expected to be relevant 
only when the shock velocity is less than about 3000 km/s, because when the relative ve- 
locity between a neutral and a proton (which immediately downstream is of order of Vgh) is 
larger than this critical value ionization becomes favored with respect to charge-exchange. 
Therefore, while during the Sedov-Taylor stage the neutral return flux has to be reckoned 
with both in working out SNR hydrodynamics and CR acceleration, it is unlikely for such a 
phenomenon to produce steep spectra during the ejecta-dominated stage as well. 

A further comment is due about the choice of the environmental parameters adopted 
in the paper. It is rather clear that the Alfvenic Mach number Ma^ oc Bq/^/pq, which 
plays a crucial role in the matter at debate, may differ quite a bit according to the nature of 
the circumstellar medium the shock propagates into. In superbubble-like environments, for 
instance, the particle density may be as low as ~ 10~^cm~^, therefore the effects discussed 
here may be even more marked. On the other hand, one could argue that for a type-la SN 
the shock should propagate in a denser and colder medium with density of order 0.1-lcm~^, 
in turn implying significantly larger Ma- 

These arguments, however, must be taken with a grain of salt, since the dependence 
on the plasma density may easily be overbalanced by a more or less efficient magnetic field 
amplification, which in many cases is expected to be mildly dependent on Bq. The absolute 
value of the downstream magnetic field, then, may be rather dependent on the actual shock 
velocity and on the circumstellar density, as discussed in commenting figure 3, therefore the 
analysis of a given remnant may differ quite a bit from this general outline. In the case of 
Tycho's SNR, for instance, assuming a particle density ~ 0.3cm~^ and = 3.7 does not 
lead to results qualitatively much different from the ones worked out here: also in that case, 
in fact, the inferred magnetic fields are strong enough to produce an appreciable steepening 
of the proton spectrum [10]. 

If any, the actual value of M^fi may be very important in the regions where magnetic 
field amplification may not be maximally efficient (like far upstream), therefore regulating 
the steepness of the spectrum up to the highest energies. In this respect, it is easy to show 
that for hot, rarefied cavities in bubble-like environments (Tq ~ 10^-10'', no ~ lO^'^cm"^) 
with standard Galactic fields Bq ~ 5//G both the sonic and the Alfvenic Mach numbers 
become of order 5-10 well within the Sedov stage. 

4 Conclusions 

In this paper we tackled the problem of particle acceleration at SNR shocks in order to 
provide a theoretical explanation for the many evidences of spectra steeper than E""^ coming 
from 7-ray observations of SNRs and from our current understanding of CR propagation in 
the Galaxy. 
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The present investigation is motivated by the fact that the prediction of the most natural 
NLDSA theory, namely that the larger the CR acceleration efficiency, the flatter the CR 
spectrum, needs to be revised in order to be consistent with current observations. 

We demonstrated that the magnetic field amplification naturally induced by the super- 
Alfvenic streaming of accelerated particles may significantly alter the properties of the CR 
scattering, actually decoupling the compression ratios felt by the fluid and by the diffusing 
particles. Under reasonable assumptions about the development and the saturation of the 
plasma instabilities (which still need to be checked against first-principle simulations given 
the impossibility of carrying out an analytical treatment of their non-linear regime), we find 
that the self-generated magnetic turbulence can lead to a steepening in the spectrum of the 
accelerated particles. 

More precisely, the expected levels of magnetic field amplification granted by streaming 
instabilities, in addition to be consistent with the fields inferred in the downstream of young 
SNRs (figure 3), can lead to CR spectra as steep as ~ E~'^-^ even in the early stages of the 
SNR evolution, i.e., when the sonic Mach number is still much larger than one. 

The effective steepening of the CR spectra turns out to be function of the CR accel- 
eration efficiency, which we tune by regulating the fraction of particles extracted from the 
thermal bath and injected in the acceleration process, rj. Such a dependence is however 
radically different from the one predicted by a NLDSA theory in which the finite velocity of 
the scattering centers is not taken into account: the spectra of the accelerated particles, in 
fact, is showed to be consistent with the test-particle prediction for very low efficiencies, but 
invariably steeper and steeper than E~'^ when rj becomes larger and larger. 

Very interestingly, in this non-linear system a larger r] on one hand produces a larger 
Per, but on the other hand it also produces a larger self-generated magnetic field, which acts 
in such a way to reduce the shock modification by steepening the CR spectrum. The net 
effect, depicted in figure 5, is that for rj > 10~^ the pressure in CRs saturates around 10-30% 
of the shock bulk pressure; also the value of the self-generated magnetic field saturates in a 
similar fashion (figure 3). This self-regulating interplay between efficient CR acceleration and 
effective magnetic field amplification may represent a key ingredient in order to quantitatively 
explain both the levels of magnetization inferred in SNRs and the CR acceleration efficiency 
required for SNRs to be the sources of Galactic CRs. 

Another important result is that the magnetic fields inferred are also the ones required to 
accelerate particles up to energies comparable with the observed knee in the diffuse spectrum 
of Galactic CRs, namely about 3-5xlO^GeV for protons (top panel of figure 4). 

We also investigated the temporal evolution of the CR spectra following the SNR during 
its ejecta-dominated and adiabatic stages, assessing the differences between instantaneous 
and cumulative spectra of the particles advected downstream (which undergo adiabatic losses 
due to the shell expansion), and of the particles escaping the system from the upstream as 
a consequence of the decrease of the SNR confining power with time [50]. We find that the 
CR acceleration efficiency is expected to drop with time, basically because of the slowing 
down of the shock due to the inertia of the swept-up material, with the net result that in 
the intermediate/late Sedov phase the SNR content in CRs is invariably dominated by the 
contributions of the earlier stages. 

This is particularly important for two reasons. First, it shows that, in order for the 
total CR contribution during the SNR lifetime to be steeper than E~'^, particles must be 
accelerated with steep spectra also during the early stages. Second, when calculating the 
expected non-thermal emission from middle-age SNRs, a time-dependent study of the SNR 
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evolution has to be carried out, since a simple snapshot of what is going on at the shock 
would easily lead to an overestimate of the spectral slope and to a underestimate of the SNR 
content in non-thermal particles. 
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